Preparation and relaxation of very stable glassy states of a simulated liquid 
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We prepare metastable glassy states in a model glass-former made of Lennard- Jones particles 
by sampling biased ensembles of trajectories with low dynamical activity. These trajectories form 
an inactive dynamical phase whose 'fast' vibrational degrees of freedom are maintained at ther- 
mal equilibrium by contact with a heat bath, while the 'slow' structural degrees of freedom are 
located in deep valleys of the energy landscape. We examine the relaxation to equilibrium and the 
vibrational properties of these metastable states. The glassy states we prepare by our trajectory 
sampling method are very stable to thermal fluctuations and also more mechanically rigid than 
low-temperature equilibrated configurations. 



As a supercooled liquid is cooled towards its glass tran- 
sition, its viscosity increases dramatically while its struc- 
ture changes only subtly [IHS]. Thus, different fluid states 
with similar structures may have relaxation times that 
differ by many orders of magnitude. In this report, we 
focus on fluid configurations that relax especially slowly. 
We do so with a field s that suppresses trajectories with 
appreciable particle motion [6HTT]. It is this field that 
controls a dynamical or space-time phase transition [7]- 
[9] in glass forming liquids, a transition between active 
fluid states and inactive states where structural relax- 
ation may be completely arrested. 

We consider a binary mixture of spherical particles 
which supports both active and inactive states. The 
structure of the inactive state differs subtly from the 
active one, and these differences render the inactive 
state extraordinarily stable. Thus, while the field s 
biases the dynamics of the system, the fluid responds 
by changing its structure, so as to arrive in long-lived 
metastable states. We find that these states are lo- 
cated in (or near [T^]) deep valleys of the energy land- 
scape [T3]. The relationships between long-lived 
metastable states and glassy behaviour have been dis- 
cussed extensively [H H4HT8] , However, even the defi- 
nition of a metastable state requires a dynamical con- 
struction that accounts for its lifetime [T51 |TS] , while the 
energy landscape is a purely static object. Since the field 
s couples directly to the dynamical evolution of the sys- 
tem, we find that it is a powerful new tool for analysing 
long-lived metastable states. 

The model we study is the Lennard- Jones (LJ) mixture 
of Kob and Andersen (KA) [19]. There are N particles 
in the system, of which Na = 0-8N are of type A and 
N B = 0.2iV are of type B. The unit of length is the 
diameter a of the type A particles, and we set the LJ 
energy for AA interactions to be e = 1. All particles 
have mass m and we take Boltzmann's constant &b = 1. 
To facilitate sampling of the s-ensemble, we consider a 
small system of N = 150 particles in a box of size (5a) 3 
with periodic boundaries, as in [pj. 



The system is coupled to a heat bath so its dynami- 
cal evolution is stochastic. We consider both Newtonian 
dynamics coupled to a thermostat, and a Monte Carlo 
(MC) dynamical scheme. Both methods give similar re- 
sults, both at equilibrium [5D] and in the s-ensemble [9J. 
We use x — (ri, r 2 , . . . , rjy) to represent the positions 
of all particles in the system. We consider ensembles of 
trajectories ('s-ensembles') based on large deviations [7 
of the dynamical activity. Within the s-ensemble, trajec- 
tories have length t \> s and the probability of a trajectory 
x(t) is 

-sK[x(t)] 

Prob[x(t)\s] = Prob[x(t)10] , (1) 

Z(s) 

where Prob[a;(t)|0] is the probability of the trajectory 
x(t) at equilibrium and Z(s) is a normalisation factor. 
The dynamical activity K measures the amount of mo- 
tion that takes place in a trajectory, and is defined by 
K = MY^AYlfjo \ri(tj + At) - r,(^)| 2 where the 
tj = jAt are equally spaced times along the trajectory, 
M = £ b s / At, and the index i runs over all particles of 
type A. The method exploits the idea that since the most 
striking glassy properties are dynamical in nature [U [5] , 
the dynamical activity is a natural order parameter for 
the glass transition [3T]. We sampled these ensembles 
using transition path sampling [5J [52] . 

We focus on inactive configurations taken from the in- 
active state in the s-ensemble, and we compare them with 
thermally-equilibrated configurations. To assess the sta- 
bility of different configurations, we used them as ini- 
tial conditions for simulations with MC dynamics, im- 
plemented as in [9j 120] . All simulations are run at tem- 
perature T = 0.6, and no biasing field s was applied. 
Results are shown in Fig. [T] where we show the mean 
square displacement of the type A particles, (r 2 (t)), and 
also their self- intermediate scattering function, F s (k,t). 
We use these simulations to model the 'melting' of the in- 
active state, and we compare this process with the heat- 
ing of a supercooled liquid state from one temperature 
to another (see also the recent experiments in |23)h In 
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FIG. 1: Self-intermediate scattering function, F s (k, t) = 
N^ 1 (Xli^i ex P (— ife • [ r i(~t) — r i(fi)])), and mean-squared dis- 
placement, {r 2 (t)) y from simulations at T — 0.6. We 
show time-dependent expectation values evaluated with equi- 
librated initial conditions at T = 0.6 (dot-dashed); from 
an inactive s-ensemble at T — 0.6 (full line, see the main 
text for details); and from equilibrated initial conditions at 
Tmit = 0.47 (dashed line). In the definition of F a (k,t), the 
sum runs over all particles of type A and k = |fc| = 7.251/ct 
corresponds to the first peak of the structure factor. 
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FIG. 2: Average energy (E), average inherent structure en- 
ergy (Eis) and average vibrational energy (Evib), for equilib- 
rium states at various temperatures, and for inactive config- 
urations. Error bars show the size of sample-to-sample fluc- 
tuations for these small systems; numerical uncertainties are 
much smaller than these error bars. In (c), the solid line is the 
result for harmonic vibrations, (i?vib) = ^NksT. (The en- 
semble of inactive configurations is the same as that in Fig. [I]) 



our MC simulations the unit of time is At, defined such 
that the diffusion constant in the limit of low density is 
Do = a 2 /At [9]. For simulations with Newtonian dy- 
namics, we take At = 1.92y 'ma 2 / 'e which allows quan- 
titative comparison with MC results. Inactive configura- 
tions were obtained from the mid-point (t — t \, s /2) of 
trajectories x(t), taken from an s-ensemble with MC dy- 
namics at T = 0.6, t obs = 150Ai and s = 0.0725/(cr 2 Ai). 
This s-ensemble is in the inactive state: we have consid- 
ered other ensembles from this state and their behaviour 
is qualitatively similar. 

For simulations with inactive initial configurations, 
(r 2 (t)) shows a plateau, with the system remaining stable 
for at least 50At before the particles diffuse away from 
their initial positions. We conclude that the inactive con- 
figurations are localized in metastable states, and must 
overcome significant free energy barriers before they re- 
lax to equilibrium. Comparing initial conditions from 
the inactive phase with equilibrated fluid configurations 
from T = 0.47, we see that these fluid states are less sta- 
ble, and relax more quickly to equilibrium. While steady 
state simulations at equilibrium and in the s-ensemble 
are similar for both MC and Newtonian dynamics, melt- 
ing and heating processes do depend significantly on the 
dynamics used in our simulations. MC dynamics approx- 
imate the overdamped limit of strong coupling to a heat 
bath, and are convenient for demonstrating the metasta- 
bility of the inactive phase, as in Fig. [T] 

In Fig.^a), we show the average energies (E) for equi- 
librated states at various temperatures, and for the in- 



active configurations. The energy of the inactive state 
is lower than the equilibrated state at the same tem- 
perature, but this difference is small compared to the 
variation in energy between different equilibrated states. 
Given that the inactive configurations are much more sta- 
ble than the thermally-equilibrated ones, their relatively 
large energy may seem surprising. 

To understand this result, we consider inherent struc- 
tures (ISs) [26], obtained by using a conjugate gradi- 
ent method to find the 'nearest' energy minimum to 
any configuration. The energy of configuration x is 
E{x) = Eis(x)+E v ih(x) where Ejs(x) is the energy of the 
inherent structure associated with x and we loosely iden- 
tify Evih(x) with 'vibrations' around the IS. Fig. [2jb,c) 
shows the averages of Eis and £vib- The inactive config- 
urations have IS energies that are lower than any of the 
equilibrated systems we considered. In computer simu- 
lations, the KA mixture has been equilibrated at tem- 
peratures as low as T — 0.42 [28]. The average inherent 
structure energy in the inactive state appears to be con- 
sistent with that of equilibrated states near to or below 
this temperature. Making the simple approximation of 
thermally-equilibrated harmonic vibrations about the IS 
positions, we predict (-Evib) = ^Nk^T, consistent with 
the data for both both thermally-equilibrated and inac- 
tive states [see Fig. |2jc)]. 

Thus, we attribute the stability of the inactive con- 
figurations (Fig. [I]) to their low inherent structure ener- 
gies. This link is consistent with studies of the energy 
landscape at equilibrium, although there is also evidence 
that slow particle motion is correlated not just with deep 
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FIG. 3: Vibrational density of states D(uj) (scaled by ui 2 ) 
for equilibrium states at T — 0.6 (dot-dashed) and T — 0.5 
(dashed), and inactive states at T = 0.6 (full line). Note 
the relative absence of low frequency modes in the inac- 
tive state. The inactive data are taken from an s-ensemble 
with Newtonian dynamics and t b s = 600A£ sampled at 
s — 0.009/(<r 2 A£), near to space-time phase coexistence, but 
restricted to K/(Nt h s a 2 ) < 0.03 [9]. Configurations were 
taken from all times throughout these trajectories. This s- 
ensemble was chosen to optimise statistics for D(uj): results 
for the inactive configurations considered in Fig.[l]are similar. 
The inset shows the participation ratio L(u). 

minima but also with saddles that have few unstable di- 
rections |3l [13l [27] . Comparing active (equilibrated) and 
inactive configurations at T = 0.6, we see from Fig. |2jb) 
that the biasing field s has a strong effect on the IS de- 
grees of freedom, while the vibrational degrees of freedom 
remain close to equilibrium at temperature T. Thus, for 
the relatively small value of s that we are considering, it 
appears that the probability of finding a configuration x 
in the inactive s-ensemble is approximately 

P(x\s) oc7>(x IS | S )e- £ « b ^/ T , (2) 

where V(xis\s) is an s-dependent statistical weight asso- 
ciated with the inherent structure Xis, while the Boltz- 
mann factor on the right hand side indicates that the 
vibrational degrees of freedom are close to equilibrium 
at the bath temperature. At equilibrium, one has 
V(xis\0) = e- Els ^ T but Fig.^b) shows that V{xi S \s) 
is dominated by ISs that are much lower in energy than 
those found at equilibrium. 

We have calculated the vibrational densities of states 
for these states by expanding the energy E(x) around 
the IS and diagonalizing the Hessian matrix to obtain 
(dimensionless) eigenfrequencies w a and eigenvectors e a . 
The density of states D(w) is the distribution of eigen- 
frequencies: eigenvectors with small u) are 'soft direc- 
tions' on the energy landscape, which may be corre- 
lated with the motion of particles during structural relax- 
ation [2H [25] . Fig. [3] shows that inactive configurations 
have fewer soft directions than configurations from ther- 
mal equilibrium: in this sense, the inactive state is more 
rigid than the thermally equilibrated states. 
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FIG. 4: Time-dependent energy in 'melting' simulations at 
T — 0.6. For low-temperature equilibrated initial conditions, 
energy flows into the system in two stages, corresponding to 
fast (t < 0.1 At) and slow (t > At) relaxation. For inactive 
initial conditions, there is only a single stage. The solid black 
line is an exponential fit with characteristic time 290At. 

We also show the participation ratio [55,, L(u) = 
(l/[NJ2i( e a ' e a) 2 ])j where the sum runs over all parti- 
cles, the vector e % a contains the components of e a associ- 
ated with particle i, and the average is over modes with 
frequency co a = u>, from all relevant configurations. In 
all cases, L(uS) decreases for small w, indicating that the 
soft modes are localized on a relatively small number of 
particles. Thus, while the inactive states have fewer soft 
directions and hence smaller vibrational fluctuations, the 
nature of the modes themselves appears similar between 
active and inactive states. 

In Fig.|4j we show the time evolution of the energy for 
the 'melting' simulations discussed above (recall Pig. fTJ . 
On taking an equilibrated configuration from T = 0.47 
and running MC dynamics at temperature T — 0.6, en- 
ergy flows into the system in two distinct stages: the 
vibrational degrees of freedom respond quickly to the 
change in temperature while the structural degrees of 
freedom respond more slowly. On the other hand, on tak- 
ing an inactive configuration and running MC dynamics 
at T = 0.6, the fast degrees of freedom in the inactive 
state are already close to equilibrium and there is no ini- 
tial stage of relaxation. The system remains localised in 
the metastable inactive state until it finally relaxes back 
to equilibrium, with an approximately exponential time- 
dependence. 

It is natural to ask what structural features of the in- 
active configurations are responsible for their low IS en- 
ergies. As in [5], we exclude crystalline states from the 
s-cnscmbles we consider, since we are specifically inter- 
ested in amorphous glassy states. Performing a common 
neighbour analysis [301 EI] j we find that inherent struc- 
tures from the inactive state are slightly richer in the 
'155' environment than their equilibrated counterparts. 
The 155 environment is associated with icosahedral co- 
ordination [31j . However, the differences are subtle and 
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sample-to-sample fluctuations large [3]: we did not find 
a specific structural motif to which we can attribute the 
stability of the inactive configurations. 

We end with a general discussion of the role of 
metastable states in the s-ensemble (see also 10J). On 
taking an initial configuration from a metastable state 
a and simulating equilibrium dynamics, the probability 
that the system remains in state a throughout a long 
time i bs is P(at — > a) ~ e -7 * obs , where 7 is a rate for re- 
laxation to equilibrium. For metastable states with long 
lifetimes, one expects a nucleation mechanism for relax- 
ation: nucleation may take place at any position in a 
large system so that 7 oc TV on taking the thermody- 
namic limit N — ¥ 00. Thus, for large enough N and i bs, 
one expects P(a — > a)/P(a — > eq) ~ e _7oAr * obs where 
P(a — > cq) w 1 is the probability of the system relaxing 
back to equilibrium. 

Let the mean dynamical activity K for long trajecto- 
ries localised in state a be k a Nt ^ s , and the mean activ- 
ity for trajectories that relax to equilibrium be k e qNt \, s . 
Then, in the s-ensemble, Eq. ([lj yields the ratio of prob- 
abilities for remaining localised in state a and for relax- 
ation to equilibrium, 



P s (a 
P s (a- 



<■ a) 
eq) 



, [s (fc cq - k a ) -70 ] Nt obB 



(3) 



where we assumed that k a and fc eq depend only weakly 
on s for small s, consistent with our observation that 
fast (vibrational, intra-state) degrees of freedom are af- 
fected weakly by s. Eq. ([3| shows that if state a is 
less active than the equilibrium state (k a < k cq ) and if 
s > s* — 7o/(fc C q~ k a ), then trajectories starting in state 
a will remain localised in that state, and will not relax to 
equilibrium even as i b s — ► 00. This construction shows 
how metastable states that are irrelevant at equilibrium 
may dominate the s-ensemble defined in ([I]). [The prob- 
ability of relaxation to a new metastable state a! ^ a 
might be larger than P(a — > a) but that is not relevant 
for the current argument.] 

The field s* required to stabilise state a may be very 
small if the metastable state is long-lived (70 is small). 
However, for small enough s, there is always a regime 
s < s* , where relaxation to equilibrium is preferred to lo- 
calisation in a metastable state, as long as 70, k a and fc oq 
are strictly positive (non-zero) constants. The definition 
of K considered here ensures that k a and k cq are both fi- 
nite. Assuming finite short-ranged interaction potentials 
and that the equilibrium state of the system is indeed a 
fluid, the nucleation rate 70 must also be non-zero even in 
the thermodynamic limit [IS]. Thus, for these systems, 
we expect any transitions in the s-ensemble to take place 
at s = s* , with s* strictly greater than zero. There are 
exceptions to this rule in idealised model systems: for 
example, in mean-field models it may be that 70 —> as 
N — > 00 due to diverging free energy barriers jlOj, while 
"kinetic constraints" can lead to 70 — ► in the thermo- 



dynamic limit [5J. Transitions at s* = might also be 
possible if the difference in activity density k cq — k a were 
to diverge, which may be relevant for glass formers |32) . 

Finally, we note that in any system with long-lived 
metastable states, the "mean-field" analysis leading to 
Eq. ([3]) predicts a dynamic phase transition at s = s*. 
However, fluctuations may destroy these transitions. For 
example, as well as P(a — > a) and P(a — ¥ eq), one 
should consider the possibility that one part of a trajec- 
tory remains localised in state a while another part has a 
structure compatible with thermal equilibrium. If this is 
likely, increasing s may result in a smooth crossover from 
active to inactive behaviour, with no dynamical phase 
transition. As demonstrated in for a kinetically con- 
strained model, it is the strength of the coupling between 
the dynamics in different parts of a system that deter- 
mines whether a dynamical phase transition takes place. 

We conclude that the s-ensemble provides a most effec- 
tive method for sampling metastable states in glassy sys- 
tems. By biasing trajectories according to their dynami- 
cal activity, the method samples these states "democrat- 
ically" , without any assumptions about their structural 
features or long-ranged correlations. In the KA mixture, 
we find metastable states that are associated with deep 
minima of the energy landscape and have few soft vi- 
brational modes. Now that these states can be prepared 
and characterised precisely, it will be interesting to see 
whether their properties can be predicted and explained 
by theories of the glass transition. 
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